
matern.estK=function (X, startpar = c(sigma2 = 1, alpha = 1), lambda = NULL, nu=1/2, 
		q = 1/4, p = 2, rmin = NULL, rmax = NULL, ...) 
{
	dataname <- deparse(substitute(X))
	if (inherits(X, "fv")) {
		K <- X
		#if (!(attr(K, "fname") %in% c("K", "Kinhom"))) 
		#	warning("Argument X does not appear to be a K-function")
	}
	#else if (inherits(X, "ppp")) {
	#	K <- Kest(X)
	#	dataname <- paste("Kest(", dataname, ")", sep = "")
	#	if (is.null(lambda)) 
	#		lambda <- summary(X)$intensity
	#}
	else stop("Unrecognised format for argument X")
	startpar <- check.named.vector(startpar, c("sigma2", "alpha"))
	Integrand <- function(r, par,nu) {
		#here is the modification
		if(nu==Inf){
			return(2*pi*r*exp(par[1]*exp(-1/2*(r/par[2])^2)))
		}#else if (nu ==1/2){
		#	return(2*pi*r*exp(par[1]*exp(-r/par[2])))
		#}
		else{
			return(2*pi*r*exp(Covariance(r,model="matern",param=c(0.0,par[1],0.0,par[2],nu))))
		}
	}
	theoret <- function(par, rvals, ..., integrand,nu) {
		if (any(par <= 0) | par[2]>5000) 
			return(rep(Inf, length(rvals)))
		th <- numeric(length(rvals))
		th[1] <- if (rvals[1] == 0) 
					0
				else integrate(integrand, lower = 0, upper = rvals[1], 
							par = par,nu=nu)$value
		for (i in 2:length(rvals)) th[i] = th[i - 1] + integrate(integrand, 
					lower = rvals[i - 1], upper = rvals[i], par = par,nu=nu)$value
		return(th)
	}
	#result <- minicontrast.weighted(K, theoret, startpar, ctrl = list(q = q, 
	#				p = p, rmin = rmin, rmax = rmax), fvlab = list(label = "%s[fit](r)", 
	#				desc = "minimum contrast fit of LGCP"), explain = list(dataname = dataname, 
	#				fname = attr(K, "fname"), modelname = "Cox process with matern pair correlation function"), 
	#		..., integrand = Integrand,nu=nu, low.weight=TRUE)
	result <-  mincontrast(K, theoret, startpar, ctrl = list(q = q, 
					p = p, rmin = rmin, rmax = rmax), fvlab = list(label = "%s[fit](r)", 
					desc = "minimum contrast fit of LGCP"), explain = list(dataname = dataname, 
					fname = attr(K, "fname"), modelname = "Cox process with matern pair correlation function"), 
			..., integrand = Integrand,nu=nu)
	
	par <- result$par
	names(par) <- c("sigma2", "alpha")
	result$par <- par
	mu <- if (is.numeric(lambda) && length(lambda) == 1 && lambda > 
					0) 
				log(lambda) - par[["sigma2"]]/2
			else NA
	result$modelpar <- c(sigma2 = par[["sigma2"]], alpha = par[["alpha"]], 
			mu = mu)
	result$internal <- list(model = "matern")
    attr(result,"nu")=nu
	return(result)
}


best.matern.estK=function(X, startpar = c(sigma2 = 1, alpha = 1), lambda = NULL, nu=c(0.25), 
		q = 1/4, p = 2, rmin = NULL, rmax = NULL, ...){
	bestmodel=list()
	minivalues=numeric()
	for(i in 1:length(nu)){
		bestmodel[[i]]=matern.estK(X, startpar,rmax=rmax,rmin=rmin,nu=nu[i])
		minivalues[i]=bestmodel[[i]]$opt$value
	}
	minione=which(minivalues==min(minivalues))[1]
	return(bestmodel[[minione]])
}